function DecodingProbability_v25(q,g)

% Clear commandline
clc
close all

% Time it
tic

% Transmitter range
T_min=30;
T_step=1;
T_max=140;

% Field cardinality is the q vector input
% Layer selection probability is the g vector input

% Layer dimensions (these are dejan examples)
k1=33;
k2=67;
K1=k1;
K2=k1+k2;

% Calculate decoding probabilities #set 1
l1_prob_1=layer_1_decod_prob(T_min,T_step,T_max,K1,K2,q(1),g);
l2_prob_1=layer_2_decod_prob(T_min,T_step,T_max,K1,K2,q(1),g);

% Calculate decoding probabilities #set 2
l1_prob_2=layer_1_decod_prob(T_min,T_step,T_max,K1,K2,q(2),g);
l2_prob_2=layer_2_decod_prob(T_min,T_step,T_max,K1,K2,q(2),g);

% Plot nice graph!
hold all

% Linespec options here
% http://www.mathworks.se/help/techdoc/ref/linespec.html

plotter(l1_prob_1,l2_prob_1,T_min,T_step,T_max,K1,K2,'r',':','o');
legend_h=plotter(l1_prob_2,l2_prob_2,T_min,T_step,T_max,K1,K2,'b',':','s');

% Set proper legend
label_1=strcat('L1, \Gamma=',num2str(g(1)),' , ',num2str(g(2)),', FF(2^1)');
label_2=strcat('L2, \Gamma=',num2str(g(1)),' , ',num2str(g(2)),', FF(2^1)');
label_3=strcat('L1, \Gamma=',num2str(g(1)),' , ',num2str(g(2)),', FF(2^8)');
label_4=strcat('L2, \Gamma=',num2str(g(1)),' , ',num2str(g(2)),', FF(2^8)');
set(legend_h,'String',[label_1;label_2;label_3;label_4],'location','SouthEast')
set(legend_h,'interpreter','tex')

% Save plot
figname=strcat('uep_ew_analytic','_g1_',num2str(g(1)),'_g2_',num2str(g(2)),'.eps');
print(gcf,'-depsc2',figname)

% Time it
toc

end

% Calculate layer 1 probability
function l1_prob = layer_1_decod_prob(T_min,T_step,T_max,K1,K2,q,g)

l1_prob=zeros(T_max,1); % Decoding 1. Layer probabilities
sol=zeros(T_max+1,1); % Array for temporary results

for tx =T_min:T_step:T_max % For each number of recv packets
    
    for n = 0:tx % For each permutation of a number of recv packets
        
        % Layer 1 by itself
        val=PM(n,K1,K1,q);
        
        val_test=0;
        
        % Layer 1 and Layer 2 gives rank K2 (This way we also get L1!)
        % Sum prob for all possible ways to achieve rank K2 with given permutation of recv packets
        
        for i=0:K1-1
             tmp1=PM(n,K1,i,q);
            
            %Optimization We only need to calculate tmp2 if tmp1!=0
            if tmp1~=0
                tmp2=PM(tx-n,K2-i,K2-i,q);
                val_test=val_test+tmp1*tmp2;
            end
            
            
        end
        
        % We have counted all the ways L1 can become full rank by itself
        % We have counted all the ways L2 can become full rank (except when
        % L1 is full rank!)
        % Both outcomes are valid for getting L1 and since they are disjoint we can just add them!
        % P(a)+P(b)-P(ab), where P(ab)=0 because they are disjoint!
        
        sol(n+1)=(val+val_test)*binopdf(n,tx,g(1));
        
    end
    
    l1_prob(tx)=sum(sol);
    sol=zeros(T_max,1);
    
    disp(['Layer 1: ' num2str(tx) ' out of ' num2str(T_max)])
    
end

end

% Calculate layer 2 probability
function l2_prob = layer_2_decod_prob(T_min,T_step,T_max,K1,K2,q,g)

% Decoding 2. Layer probabilities
l2_prob=zeros(T_max,1);

% Array for temporary results
sol=zeros(T_max+1,1);

for tx =T_min:T_step:T_max % For each number of recv packets
        
    for n = 0:tx % For all permutations of recv packets
        
        val=0;
        
        for i=0:K1 % For all possible ways to achieve rank K2 with given permutation of recv packets
            
            tmp1=PM(n,K1,i,q);
            
            % Optimization no need for tmp2 when tmp1=0!
            if tmp1==0
                continue;
            end
            
            tmp2=PM(tx-n,K2-i,K2-i,q);
            val=val+tmp1*tmp2;
        end
        
        sol(n+1)=val*binopdf(n,tx,g(1));
        
    end
    
    l2_prob(tx)=sum(sol);
    sol=zeros(T_max,1);
    
    disp(['Layer 2: ' num2str(tx) ' out of ' num2str(T_max)])
    
end

end

% Plotter for a nice graph!
function legend_h = plotter(l1_prob,l2_prob,T_min,T_step,T_max,K1,K2,color,style,symbol)

% Replace 0 with NaN in (l1_prob,l2_prob) for prettier plot
for k=1:length(l1_prob)
    
    if l1_prob(k)==0
        l1_prob(k)=NaN;
    end
end
for k=1:length(l2_prob)
    if l2_prob(k)==0
        l2_prob(k)=NaN;
    end
    
end

% Plotting
figure(1)
plot(1:length(l1_prob),l1_prob,'-*','Color',color,'LineStyle','-','Marker',symbol)
plot(1:length(l2_prob),l2_prob,'-*','Color',color,'LineStyle',style,'Marker',symbol)

% Plot annotation
grid('on')
% pbaspect([2.5 1 1])
legend_h = legend('location','SouthEast');
set(gca,'XTick',0:10:T_max)
set(gca,'YTick',0:0.1:1)
xlim([T_min T_max])
ylim([0 1])

end

% The New helper function
function P = PM(m,n,r,q)

P1=zeros(n+1,n+1);
for i=1:length(P1)
    entry_val=1/(q^(n-(i-1)));
    P1(i,i)=entry_val;
    if i<n+1
        P1(i+1,i)=1-entry_val;
    end 
end

s1=zeros(n+1,1);
s1(1)=1;

val=(P1^m)*s1;
P=val(r+1);

end




